(*  :Title:	Discrete Fourier Transforms  *)

(*  :Authors:	Wally McClure, Brian Evans, James McClellan  *)

(*
    :Summary:	This file allows the computation of the forward and
		inverse discrete Fourier transform (DFT).
 *)

(*  :Context:	SignalProcessing`Digital`DFT`  *)

(*  :PackageVersion:  2.7	*)

(*
    :Copyright:	Copyright 1989-1991 by Brian L. Evans
		Georgia Tech Research Corporation

	Permission to use, copy, modify, and distribute this software
	and its documentation for any purpose and without fee is
	hereby granted, provided that the above copyright notice
	appear in all copies and that both that copyright notice and
	this permission notice appear in supporting documentation,
	and that the name of the Georgia Tech Research Corporation,
	Georgia Tech, or Georgia Institute of Technology not be used
	in advertising or publicity pertaining to distribution of the
	software without specific, written prior permission.  Georgia
	Tech makes no representations about the suitability of this
	software for any purpose.  It is provided "as is" without
	express or implied warranty.
 *)

(*
    :History:   start			July 30, 1989
		redirection		October 2, 1989
		made into package	April 18, 1990
		made into release	May 25th, 1990
 *)

(*  :Keywords:	discrete Fourier transform  *)

(*
    :Source:	'Discrete-Time Signal Processing', Alan V. Oppenheim
		 	and Ronald W. Schafer
 *)

(*  :Warning:	*)

(*  :Mathematica Version:  1.2 or 2.0  *)

(*  
    :Limitation: The DFT is calculated by first attempting to calculate the 
    		 DTFT of a function using the MyDTFT rule base.  At the end 
		 of the MyDTFT rule base, a call to the z-transform is made.  
		 If  the  function  does  not have a transform specified in 
		 either domain, then a nasty expression is returned.    
 *)

(*
    :Discussion:  The DFT is based on the DTFT.  A signal is made finite
		  in extent by multiplying it by ( Step[n] - Step[n - N] ),
		  where N is the number of points at which to evaluate the
		  DFT.
 *)

(*
    :Functions:	DFTransform
		InvDFTransform
 *)



(*  B E G I N     P A C K A G E  *)

BeginPackage[ "SignalProcessing`Digital`DFT`",
	      "SignalProcessing`Digital`DTFT`",
	      "SignalProcessing`Digital`DSupport`",
	      "SignalProcessing`Digital`InvZTransform`",
	      "SignalProcessing`Digital`ZTransform`",
	      "SignalProcessing`Digital`ZSupport`",
	      "SignalProcessing`Support`TransSupport`",
	      "SignalProcessing`Support`ROC`",
	      "SignalProcessing`Support`SigProc`",
	      "SignalProcessing`Support`SupCode`",
	      If [ TrueQ[$VersionNumber >= 2.0],
		   "Algebra`SymbolicSum`",
		   "Algebra`GosperSum`" ] ]


If [ TrueQ[ $VersionNumber >= 2.0 ],
     Off[ General::spell ];
     Off[ General::spell1 ] ];


(*  U S A G E     I N F O R M A T I O N  *)

DFTData::usage =
	"Data-tag for the symbolic DFT."

DFTransform::usage =
	"DFTransform[function, NumberOfPoints, TimeVariables, \
	FourierVariables, options] returns the DFT of function, \
	where function is a function of TimeVariables defined from \
	0 to NumberOfPoints - 1. \
	Note that only the first two arguments \
	are required and that DFTransform calls DTFTransform."

Finish::usage =
	"Finish is a data-tag for the DFT."

InvDFTransform::usage =
	"InvDFTransform[function, NumberOfPoints, FourierVariables, \
	TimeVariables, options] takes the inverse DFT of function, \
	where function is a function of FourierVariables defined from \
	0 to NumberOfPoints - 1. \
	Note that only the first two arguments \
	are required and that InvDFTransform calls InvDTFTransform."

KVariables::usage =
	"KVariables is a data-tag for the DFT."

Start::usage =
	"Start is a data-tag for the DFT."

(*  E N D     U S A G E     I N F O R M A T I O N  *)


Begin[ "`Private`" ]


(*  displayFixUp  *)
displayFixUp[ trans_, op_ ] :=
	Append[ fixUp[op, trans],
		TransformLookup -> Replace[TransformLookup, op] ]

(*  validvarQ  *)
validvarQ[ x_Symbol ] := True
validvarQ[ x_List ] := Apply[And, Map[VariableQ, x]]
validvarQ[ x_ ] := False


(*  B E G I N     D F T  *)

(*  Evaluation of DFT operator  *)
Unprotect[DFT]
DFT/: TheFunction[ DFT[N_, n_, k_] [f_] ] := DFTransform[f, N, n, k]
Protect[DFT]

(*  Options for DFTransform  *)
DFTransform/: Options[ DFTransform ] := Options[ DTFTransform ]

(*  Extension of TheFunction to return the transform function.  *)
DFTData/: TheFunction[ DFTData[ trans_, var_ ] ] := trans

(*  Magnitude/Phase plot of a DFT object *)
mppoptions = {	Domain -> Discrete, DomainScale -> Linear,
		MagRangeScale -> Linear, PhaseRangeScale -> Degree,
		PlotRange -> All }

DFTData/: MagPhasePlot[ DFTData[trans_, KVariables[k_Symbol],
				Start[st_], Finish[end_] ] ] :=
	MagPhasePlot[ trans, {k, 0, end - st}, mppoptions ]

DFTData/: MagPhasePlot[ DFTData[trans_, KVariables[{k1_Symbol, k2_Symbol}],
				Start[{s1_, s2_}], Finish[{e1_, e2_}] ] ] :=
	MagPhasePlot[ trans, {k1, 0, e1 - s1}, {k2, 0, e2 - s2}, mppoptions ]

(*  This Code decides between a one-dimensional or multidimensional DFT. *)
DFTransform[ f_ ] :=
	Message[ Transform::novariables, "N (length)", GetVariables[f] ]

DFTransform[ f1_, N_ ] :=
	DFTransform[ f1, N, DummyVariables[Length[N], Global`n] ]

DFTransform[ f1_, N_, n_ ] :=
	DFTransform[ f1, N, n, DummyVariables[Length[n], Global`k] ]

DFTransform[ fun_, N_, n_, k_, options___ ]:=
	Message[ Transform::badvar, "discrete-time", DFTransform, n ] /;
	! validvarQ[n]

DFTransform[ fun_, N_, n_, k_, options___ ]:=
	Message[ Transform::badvar, "discrete-frequency", DFTransform, k ] /;
	! validvarQ[k]

DFTransform[ f1_, N_, n_, k_, options___ ] :=
	Block [	{begin, end, newfun, notenough, op, trans, vars, w, wvars},
		notenough = ( Length[k] < Length[n] );
		If [ notenough,
		     Message[Transform::notenough, "k (DFT index)"] ];
		vars = If [ notenough, DummyVariables[Length[n], Global`k], k ];
		wvars = DummyVariables[Length[n], w];

		op = ToList[options] ~Join~ Options[DFTransform];
		trans = If [ Length[n] > 1,  
			     multiDDFT[f1, N, n, vars, wvars, op],
			     MyDFT[f1, N, n, vars, wvars, op] ];

		begin = Start[ If [AtomQ[n], 0, Table[0, {Length[n]}]] ];
		end = Finish[ N - 1 ]; 

		DFTData[FourierSimplify[trans], begin, end, KVariables[vars]] ]

(*  multiDDFT --  multidimensional DFT  *)
multiDDFT[ fun_, Num_, varin_, varout_, wvars_, op_ ] :=
	Block [	{F2 = fun, jj, length, num, var, w},
		length = Length[varin];
		For [ jj = 1, jj <= length, jj++,
		      var = varin[[jj]];
		      num = Num[[jj]];
		      w = wvars[[jj]];
		      F2 = MyDFT[F2, num, var, varout[[j]], w, op] ];
		F2 ]

(*  MyDFT  *)
MyDFT[ x_, N_, n_, k_, w_, op_ ] :=
	DualDFT[ x, N, n, k, w, op, DFTransform ]

(*  E N D     D F T  *)


(*  B E G I N     I N V E R S E     D F T  *)

(*  Messages  *)
InvDFTransform::badlength = "Conflicting lengths in inverse DFT: `` != ``."

(*  Evaluation of inverse DFT operator  *)
Unprotect[InvDFT]
InvDFT/: TheFunction[ InvDFT[N_, k_, n_] [f_] ] := InvDFTransform[f, N, k, n]
Protect[InvDFT]

(*  Options for InvDFTransform  *)
InvDFTransform/: Options[ InvDFTransform ] :=
	displayFixUp[ InvDFTransform,
		      { Definition -> False } ~Join~ Options[InvDTFTransform] ]

(*  Choose whether to do the one or multi dimensional inverse DFT`s.  *)
InvDFTransform[ DFTData[x_, Start[st_], Finish[end_], KVariables[k_]] ] :=
	InvDFTransform[ x, end - st + 1, k ]

InvDFTransform[ f_ ] :=
	Message[ Transform::novariables, "N (length)", GetVariables[f] ]

InvDFTransform[ DFTData[x_, Start[st_], Finish[end_], KVariables[k_]], N_ ] :=
	If [ SameQ[ end - st + 1, N ],
	     InvDFTransform[ x, N, k ],
	     Message[ InvDFTransform::badlength, N, end - st + 1 ] ]

InvDFTransform[ f_, N_ ] :=
	InvDFTransform[ f, N, DummyVariables[Length[N], Global`k] ]

InvDFTransform[ f_, N_, k_ ] :=
	InvDFTransform[ f, N, k, DummyVariables[Length[k], Global`n] ]

InvDFTransform[ fun_, N_, k_, n_, options___ ]:=
	Message[ Transform::badvar, "discrete-frequency", InvDFTransform, k ] /;
	! validvarQ[k]

InvDFTransform[ fun_, N_, k_, n_, options___ ]:=
	Message[ Transform::badvar, "discrete-time", InvDFTransform, n ] /;
	! validvarQ[n]

InvDFTransform[ fun_, N_, k_, n_, options___ ]:=
	Block [	{notenough, nvars, op, w, wvars}, 
		notenough = ( Length[n] < Length[k] );
		If [ notenough, Message[Transform::notenough, "k (frequency)"]];
		nvars = If [ notenough,
			     DummyVariables[Length[n], Global`n],
			     n ];
		op = ToList[options] ~Join~ Options[DFTransform];
		wvars = DummyVariables[Length[n], w];

		If [ Length[k] > 1, 
		     multiDInvDFT[ fun, N, k, nvars, wvars, op ],
		     MyInvDFT[ fun, N, k, nvars, wvars, op ] ] ]

(*   Multidimensional inverse DFT  *)
multiDInvDFT[ fun_, N_, k_, n_, w_, op_ ] :=
	Block [	{F1 = fun, jj, length},

		length = Length[Num];
		For [ jj = 1, jj <= length, jj++,
		      F1 = MyInvDFT[ F1, N[[jj]], k[[jj]],
				     n[[jj]], w[[jj]], op ] ];

		F1 ]

(*  MyInvDFT  *)
MyInvDFT[ x_, N_, k_, n_, w_, op_ ] :=
	DualDFT[ x, N, k, n, w, op, InvDFTransform ]

(*  E N D     I N V E R S E     D F T  *)


(*  B E G I N     D U A L     D F T  *)

DualDFT[ x_, len_, n_, k_, w_, op_, flag_ ] :=
	Block [	{newdualrules, newexpr, newx = x, oldexpr = Null, ret, trace},

		(* generate the DFT rules *)
		newdualrules = TransformFixUp[ n, k, op, dualDFT, False,
					       DFTransform, Null, Null ];
		If [ SameQ[flag, InvDFTransform],
		     newdualrules = newdualrules /. k -> -k ];
		DualRules = Join[ newdualrules, OriginalDualRules ];

		(* determine the level of dialogue *)
		trace = SameQ[ Replace[Dialogue, op], All ];

		(* rewrite trig functions in terms of exponentials *)
		(* and put exponentials in a cannonical form       *)
		newx = x /. TrigToExpRules;
		newx = newx /. ( Exp[a_] :> Exp[ Factor[a] ] );
		If [ trace && ! SameQ[x, newx],
		     Print[x];
		     Print["which is rewritten in a more convenient exponential form as"];
		     Print[newx] ];

		(* check for inverse DFT *)
		newexpr = dualDFT[newx, n, k, w, len, flag, fixUp[op, flag]];
		If [ SameQ[flag, InvDFTransform],
		     If [ trace,
			  Print[ "Using the Duality Theorem, the inverse DFT ",
			  	 "of length ", len, " is the DFT of length ",
			  	 len, " divided by ", len, " with the index ",
			  	 "variable ", k, " negated in the result:" ] ];
		      newexpr = Rev[k][newexpr / len] ];

		(* take the 1-D dual DFT *)
		While [ ! SameQ[newexpr, oldexpr],
			If [ trace, Print[newexpr]; Print[ "which becomes" ] ];
			oldexpr = newexpr;
			newexpr = oldexpr /. ( dualDFT[a__] :>
				    Replace[dualDFT[a], DualRules] ) ];

		(* fix up the result *)
		newexpr = posttransform[oldexpr];
		While [ ! SameQ[newexpr, oldexpr],
			If [ trace, Print[newexpr]; Print[ "which becomes" ] ];
			oldexpr = newexpr;
			newexpr = posttransform[oldexpr] ];

		newexpr ]

(*  circularShiftValue  *)
n0checkvalue = Null

n0checkQ[n0_] :=
	Block [	{},
		If [ ValueQ[n0checkvalue], n0checkvalue = n0 ];
		SameQ[n0, n0checkvalue] ]

circularShiftQ[x_, n_, L_] :=
	Block [	{demodx, modtoken},
		Unset[n0checkvalue];
		demodx = x /. ( Mod[n + n0_, L] :> modtoken /; n0checkQ[n0] );
		FreeQ[demodx, n] ]

circularShiftValue[x_, n_, L_] := n0checkvalue /; circularShiftQ[x, n, L]

(*  fixUp  *)
fixUp[op_, DFTransform] :=
	{ Definition -> Replace[Definition, op],
	  Dialogue -> Replace[Dialogue, op] }

fixUp[op_, InvDFTransform] :=
	{ Definition -> Replace[Definition, op],
	  Dialogue -> Replace[Dialogue, op],
	  Terms -> Replace[Terms, op] }

(*  posttransform  *)
posttransform[x_] := x /. postrules

postrules = {

	Rev[v_][x_] :>
		If [ Count[{x}, f_[p___][a___], Infinity] > 0,
		     Rev[v][posttransform[x]],          (* operators present *)
		     posttransform[x] /. v -> -v ] ,    (* no operators *)

	Rev[k_][ Rev[k_] [ x_ ] ] :> x ,

	CircularShift[n0_, L_, n_][a_. Impulse[n_ + b_.]] :>
		a Impulse[n + b + n0] /;
		FreeQ[{a, b}, n] ,

	dualDFT[ x_, n_, k_, w_, L_, DFTransform, op_ ] :>
		Block [	{dtft},
			dtft = DTFTransform[ x Pulse[L, n], n, w, op ];
			RewriteSummations[ TheFunction[dtft], w, 0, 2 Pi ] /.
				w -> (2 Pi k / L) ],

	dualDFT[ x_, k_, n_, w_, L_, InvDFTransform, op_ ] :>
		InvDTFTransform[x /. (k -> w L / (2 Pi)), w, n, op] ,

	Impulse[-v_Symbol] :> Impulse[v]
}


Format[ dualDFT[ x_, n_, k_, w_, L_, rest___ ] ] :=
	SequenceForm[ ColumnForm[ { "DFT",
		         StringJoin["   ", ToString[L], ", ", ToString[n]] } ],
		      { x } ]

Format[ w ] := "w"


DualRules = { }

OriginalDualRules = {

(*  Constant value--- some call it noise  *)
dualDFT[ a_, n_, k_, w_, L_, flag_, op_ ] :> a L Impulse[k] /; FreeQ[a, n] ,

(*  Impulse  *)
dualDFT[ a_. Impulse[n_ + n0_.], n_, k_, w_, L_, flag_, op_ ] :>
	a Exp[-2 I Pi n0 k / L] /;
	FreeQ[{a,n0}, n] && Implies[NumberQ[n0], IntegerQ[n0]] ,

(*  Multiplication by a complex exponential  *)
dualDFT[ Exp[Complex[0, k01_] Pi k02_. (a_. n_ + b_.)] x_., n_, k_,
		w_, L_, flag_, op_ ] :>
	Block [	{cshift, dft},
		cshift = a k01 k02 L / 2;
		If [ ! IntegerQ[cshift] && SameQ[Replace[Dialogue, op], All],
		     Print[ "assuming ", cshift, " is an integer" ] ];
		dft = dualDFT[ x, n, k, w, L, flag, op ];
		Exp[-2 I Pi (b/a) k / L] CircularShift[cshift, L, k][dft] ] /;
	FreeQ[{a, b, k01, k02}, n] ,

(*  Circular shift of sequence  *)
dualDFT[ CircularShift[n0_, L_, n_][x_], n_, k_, w_, L_, flag_, op_ ] :>
	Exp[2 Pi I k n0 / L] dualDFT[ x, n, k, w, L, flag, op ] /;
	FreeQ[{n0, L}, n] ,

dualDFT[ x_, n_, k_, w_, L_, flag_, op_ ] :>
	Block [	{circShift},
		circShift = circularShiftValue[x, n, L];
		newx = x /. ( Mod[n + circShift, L] -> n );
		Exp[2 Pi I k circShift / L] *
		  dualDFT[ newx, n, k, w, L, flag, op ] ] /;
	circularShiftQ[x, n, L] ,

(*  Conjugate of sequence  *)
dualDFT[ Conjugate[x_], n_, k_, w_, L_, flag_, op_ ] :>
	Conjugate[
	  CircularShift[0, N, k][
	    Rev[k][dualDFT[x, n, k, w, L, flag, op]] ] ] ,

(*  Reverse of a sequence  *)
dualDFT[ Rev[n_][x_], n_, k_, w_, L_, flag_, op_ ] :>
	Rev[k][ dualDFT[x, n, k, w, L, flag, op] ] ,

(*  Homogeneity  *)
dualDFT[ a_ x_, n_, k_, rest__ ] :>
	a dualDFT[ x, n, k, rest ] /;
	FreeQ[a, n] ,

(*  Additivity  *)
dualDFT[ x_ + y_, n_, k_, rest__ ] :>
	dualDFT[ x, n, k, rest ] + dualDFT[ y, n, k, rest ] ,

(*  Apply the Definition  *)
dualDFT[ x_, n_, k_, w_, L_, flag_, op_ ] :>
Block [ {newx, state, sumhead, trans},
	newop = fixUp[ Prepend[op, Definition -> False] ]; (* disable option *)
	sumhead = If [ TrueQ[$VersionNumber >= 2.0], SymbolicSum, GosperSum ];
	newx = ToDiscrete [ TheFunction[x] ];
	trans = sumhead[newx Exp[-2 Pi n k / L], {n, 0, L-1}];
	If [ SameQ[Head[trans], sumhead],
	     dualDFT[x, n, k, w, L, flag, newop],
	     TransformDialogue[Definition, op, sumhead, x, trans] ] ] /;
	Replace[Definition, op]

}

(*  E N D     D U A L     D F T  *)





(*  T R I G     R E W R I T E     R U L E S  *)

(* Trigonometric to Complex Exponential Simplification rules *)

TrigToExpRules = {
	Sin[ a__ ]	:> 	( Exp[ I a ] - Exp[ - I a ] ) / ( 2 I ),
	1/Sin[ a__ ]	:> 	1/( ( Exp[ I a ] - Exp[ - I a ] ) / ( 2 I ) ) ,
	Cos[ a__ ]	:>	( Exp[ I a ] + Exp[ - I a ] ) / 2,
	1/Cos[ a__ ]	:>	1/( ( Exp[ I a ] + Exp[ - I a ] ) / 2 ),
	Tan[ a__ ]	:>	Sin[ a ] / Cos[ a ],
	Cot[ a__ ]	:>	Cos[ a ] / Sin[ a ],
	Sec[ a__ ]	:>	1 / Cos[ a ],
	Csc[ a__ ]	:>	1 / Sin[ a ]
}


(* Exponential rewrite rules *)

ExpToTrigRules = {
	(a_. Exp[ Complex[0, b_] w_. ] + a_. Exp[ Complex[0, c_] w_. ] + x_.) :>
		2 a Cos[Abs[b] w] + x /;
		( b == -c ),

	(a_. Exp[ Complex[0, b_] w_. ] - a_. Exp[ Complex[0, c_] w_. ] + x_.) :>
		2 a Sin[b w] + x /;
		( b == -c ),

	(d_ - d_. Exp[ Complex[0, b_] c_. ]) :>
		-2 I d Exp[ I b c / 2 ] Sin[ b c / 2 ],

	(d_. Exp[ Complex[0, b_] c_. ] - d_) :>
		2 I d Exp[ I b c / 2 ] Sin[ b c / 2 ],

	(d_ + d_. Exp[ Complex[0, b_] c_. ]) :>
		2 d Exp[ I b c / 2 ] Cos[ b c / 2 ]
}

(*  E N D     T R I G     R E W R I T E     R U L E S  *)


(*  E N D     P A C K A G E  *)

End[]
EndPackage[]

If [ TrueQ[ $VersionNumber >= 2.0 ],
     On[ General::spell ];
     On[ General::spell1 ] ];


(*  H E L P     I N F O R M A T I O N  *)

Combine [ SPfunctions, { DFTransform, InvDFTransform } ]
Protect [ DFTransform, InvDFTransform ]


(*  E N D I N G     M E S S A G E  *)

Print[ "The discrete Fourier transform (DFT) rule bases are loaded." ]
Null
